rm(list=ls())
gc()

library(estimatr)
library(tidyverse)
library(ggthemes)
library(cowplot)
library(scales)
library(lfe)
library(lsmeans)
library(multcompView)
library(multcomp)

# Theme for pretty plots
devtools::source_url("https://kylepeyton.github.io/assets/theme_kyle.R")

# Make some adjustments ...
theme_kyle2 <- function () { 
  theme_kyle(font_size = 22) + #%+replace%
    theme(legend.position = "bottom",
          legend.justification = "center",
          text = element_text(family = "serif", size = 18),
          plot.title = element_text(margin = margin(b = 0, l = -10), 
                                    size = 24, face = "bold", family = "serif"), 
          legend.key.size = unit(3, "line"),
          legend.margin = margin(0,0,0,0),
          legend.box.margin = margin(-30, -10, -10, -10),
          strip.text.x = element_text(size = 22, family = "serif"),
          axis.title.y = element_text(size = 19, family = "serif"),
          legend.title = element_blank(),
          axis.text.x = element_text(size = 19, family = "serif"),
          strip.background = element_rect(fill = "grey90"),
          axis.line.y = element_line(size = I(4/12), colour = "black"),
          axis.ticks = element_line(colour = "black"),
          axis.ticks.y = element_line(size = I(4/12), colour = "black"))
}

#### Read in data ----------------------

input_path <- "" # Specify where the data lives

## Create sub-directory for output:
dir.create(file.path(input_path, "output/"))


## Read in Study 1 data, subset to white subjects and black/white interactions
mturk_ug <- 
  haven::read_dta(paste0(input_path, "Study1_Deidentified.dta")) %>% 
  filter(dem_racewhite == 1, proposer_race != "other") %>% 
  mutate(sample = "MTurk",
         proposer_race = factor(proposer_race, levels = c("white", "black"))) %>% 
  haven::zap_labels()

## Read in Study 2 data, subset to whites that completed >= 75% of assigned
## rounds. 
evals_ssi <- 
  haven::read_dta(paste0(input_path, "Study2_Deidentified.dta")) %>% 
  filter(dem_racewhite == 1) %>% 
  mutate(sample = "SSI", 
         responder_race = haven::as_factor(responder_race),
         proposer_race = haven::as_factor(proposer_race),
         round_racepair = haven::as_factor(round_racepair),
         round_racepair = relevel(round_racepair, ref = "WW")) %>% 
  droplevels() 

## Eliminate people who stop answering fairness evaluations less than 75% into 
## the survey, restrict to BW interactions and eliminate folks who failed 
## attention check twice
evals_ssi <- 
  evals_ssi %>% 
  group_by(pid) %>% 
  mutate(obs_fairness = sum(!is.na(offer_fairness))) %>% 
  filter(obs_fairness != 0) %>% 
  mutate(failtwo = ifelse(is.na(failtwo), 0, failtwo),
         r = ifelse(!is.na(skip_fairness), round, NA),
         maxround = max(r, na.rm = TRUE)) %>% 
  filter(maxround > 41*.75, proposer_race != "O", responder_race != "O",
         failtwo != 1) %>% 
  ungroup() %>% 
  haven::zap_labels()

## Code up 7-point party identification scale:
mturk_ug <- 
  within(mturk_ug, {
    dem_pid7 <- rep(NA, nrow(mturk_ug))
    dem_pid7[!is.na(dem_partyid)] <- 4
    dem_pid7[dem_lean == 2] <- 3
    dem_pid7[dem_partyid == 2] <- 2
    dem_pid7[dem_strongd == 1] <- 1
    dem_pid7[dem_lean == 1] <- 5
    dem_pid7[dem_partyid == 1] <- 6
    dem_pid7[dem_strongr == 1] <- 7
  })

evals_ssi <- 
  within(evals_ssi, {
    dem_pid7 <- rep(NA, nrow(evals_ssi))
    dem_pid7[!is.na(dem_partyid)] <- 4
    dem_pid7[dem_lean == 2] <- 3
    dem_pid7[dem_partyid == 2] <- 2
    dem_pid7[dem_strongd == 1] <- 1
    dem_pid7[dem_lean == 1] <- 5
    dem_pid7[dem_partyid == 1] <- 6
    dem_pid7[dem_strongr == 1] <- 7
  })


#####-------------Study 1 descriptives (reported in manuscript)------------#####

# Number of white responders (subjects)
length(unique(mturk_ug$responder_id))

# Acceptance rate
mean(mturk_ug$accept_offer, na.rm = TRUE)

# Proportion of resentful respondents
mean(mturk_ug$symbolic_binary[!duplicated(mturk_ug$responder_id)])

# Proportion of prejudiced respondents
mean(mturk_ug$overt_binary[!duplicated(mturk_ug$responder_id)])

# Tabulate indicators for both
table(mturk_ug$overt_binary[!duplicated(mturk_ug$responder_id)],
      mturk_ug$symbolic_binary[!duplicated(mturk_ug$responder_id)])

# Classification rates by partisanship
mturk_ug %>% 
  distinct(., responder_id, .keep_all = TRUE) %>% 
  filter(dem_pid7 != 4) %>% 
  group_by(dem_pid7 >= 5) %>% 
  summarise(prejudice = mean(overt_binary, na.rm = T),
            resentful = mean(symbolic_binary, na.rm = T))

## Are Republicans more likely to discriminate than Democrats?
fit_pid <- felm(accept_offer ~ proposer_race*as.numeric(dem_pid7 >= 5) |
                  factor(offer_amount) + factor(round) | 0 | responder_id, 
                data = mturk_ug %>% filter(dem_pid7 != 4))

tidy(fit_pid) %>% mutate_if(is.numeric, round, 3)

#####-------------Study 2 descriptives (reported in manuscript)-----------#####

# Number of white responders (subjects)
length(unique(evals_ssi$pid))

# Correlation bt likelihood of acceptance and perceived fairness
cor(evals_ssi$offer_fairness, evals_ssi$offer_accepted, use = "complete.obs")

# Prejudice v. resentment proportions
evals_ssi %>% 
  distinct(., pid, .keep_all = TRUE) %>% 
  summarise(prejudice = mean(overt_binary, na.rm = TRUE),
            resentful = mean(symbolic_binary, na.rm = TRUE))

# Correlation bt measures
cor(evals_ssi$symbolic_binary, evals_ssi$overt_binary)

# Classification rates by partisanship
evals_ssi %>% 
  distinct(., pid, .keep_all = TRUE) %>% 
  filter(dem_pid7 != 4) %>% 
  group_by(dem_pid7 >= 5) %>% 
  summarise(prejudice = mean(overt_binary, na.rm = TRUE),
            resentful = mean(symbolic_binary, na.rm = TRUE))

#####-------------Additional Descriptives in Supplementary Materials-------#####

## Reliability of prejudice/resentment measures:

# Resentment, Study 1
mturk_ug %>% 
  distinct(., responder_id, .keep_all = TRUE) %>% 
  dplyr::select(dem_ks1, dem_ks2, dem_ks3, dem_ks4) %>% 
  psych::alpha()

# Prejudice, Study 1
mturk_ug %>% 
  distinct(., responder_id, .keep_all = TRUE) %>% 
  mutate(bw_lazygap = dem_lazy_black - dem_lazy_white,
         bw_unintgap = dem_unint_black - dem_unint_white,
         bw_untrstgap = dem_untrst_black - dem_untrst_white,
         bw_violgap = dem_viol_black - dem_viol_white) %>% 
  dplyr::select(starts_with("bw_")) %>% 
  psych::alpha()

# Resentment, Study 2
evals_ssi %>% 
  distinct(., pid, .keep_all = TRUE) %>% 
  dplyr::select(dem_ks1, dem_ks2, dem_ks3, dem_ks4) %>% 
  psych::alpha()

# Prejudice, Study 2
evals_ssi %>% 
  distinct(., pid, .keep_all = TRUE) %>% 
  mutate(bw_lazygap = dem_lazy_black - dem_lazy_white,
         bw_unintgap = dem_unint_black - dem_unint_white,
         bw_untrstgap = dem_untrst_black - dem_untrst_white,
         bw_violgap = dem_viol_black - dem_viol_white) %>% 
  dplyr::select(starts_with("bw_")) %>% 
  psych::alpha()

# Mean scores on resentment/prejudice scales in both studies:
mturk_ug %>% 
  distinct(., responder_id, .keep_all = TRUE) %>% 
  summarise(avg_resentment = mean(symbolic_index, na.rm = TRUE),
            avg_prejudice = mean(overt_index, na.rm = TRUE),
            prop_resent = mean(symbolic_binary, na.rm = TRUE),
            prop_prejudice = mean(overt_binary, na.rm = TRUE))

evals_ssi %>% 
  distinct(., pid, .keep_all = TRUE) %>% 
  summarise(avg_resentment = mean(symbolic_index, na.rm = TRUE),
            avg_prejudice = mean(overt_index, na.rm = TRUE),
            prop_resent = mean(symbolic_binary, na.rm = TRUE),
            prop_prejudice = mean(overt_binary, na.rm = TRUE))


# Correlation between resentment/prejudice scales in both studies:
cor(mturk_ug$symbolic_index[!duplicated(mturk_ug$responder_id)], 
    mturk_ug$overt_index[!duplicated(mturk_ug$responder_id)])

cor(evals_ssi$symbolic_index[!duplicated(evals_ssi$pid)], 
    evals_ssi$overt_index[!duplicated(evals_ssi$pid)])

# Correlation w/ partisanship in stacked dataset
stacked_surveys <- 
  evals_ssi %>% 
  distinct(pid, .keep_all = TRUE) %>% 
  mutate(Study = "Study 2") %>% 
  bind_rows(mturk_ug %>% 
              distinct(responder_id, .keep_all = TRUE) %>% 
              mutate(Study = "Study 1"))

cor(stacked_surveys$dem_pid7, stacked_surveys$symbolic_index, 
    use = "complete.obs")

cor(stacked_surveys$dem_pid7, stacked_surveys$overt_index, 
    use = "complete.obs")

# Correlation w/ voting Republican in 2012 election
cor(as.numeric(stacked_surveys$dem_obama == 2), 
    stacked_surveys$symbolic_index, use = "complete.obs")

cor(as.numeric(stacked_surveys$dem_obama == 2), 
    stacked_surveys$overt_index, use = "complete.obs")

#####-------------Study 1 Fig. 1 -----------------------------------------#####
fit1_ug <- felm(accept_offer ~ proposer_race | factor(offer_amount)  + 
                  factor(round) | 0 | responder_id, data = mturk_ug)


fit2_ug <- felm(accept_offer ~ proposer_race*symbolic_binary |
                  factor(offer_amount) + factor(round) | 0 | 
                  responder_id, data = mturk_ug)

fit3_ug <- felm(accept_offer ~ proposer_race*overt_binary | 
                  factor(offer_amount) + factor(round) | 0 | 
                  responder_id, data = mturk_ug)

fit4_ug <- felm(accept_offer ~ proposer_race*overt_binary + 
                  proposer_race*symbolic_binary | factor(offer_amount) + 
                  factor(round) | 0 | responder_id, data = mturk_ug)

est_df <- 
  bind_rows(tidy(fit1_ug) %>% mutate(model = "Model 1"),
            tidy(fit2_ug) %>% mutate(model = "Model 2"),
            tidy(fit3_ug) %>% mutate(model = "Model 3"),
            tidy(fit4_ug) %>% mutate(model = "Model 4")) %>% 
  filter(!(term %in% c("overt_binary", "symbolic_binary"))) %>% 
  mutate(term = factor(term, levels = c("proposer_raceblack",
                                        "proposer_raceblack:symbolic_binary",
                                        "proposer_raceblack:overt_binary"),
                       labels = c("Black Proposer", 
                                  "Black Proposer x RRI",
                                  "Black Proposer x EPI"))) %>% 
  mutate(li90 = estimate - 1.645*std.error,
         ui90 = estimate + 1.645*std.error,
         li95 = estimate - 1.96*std.error,
         ui95 = estimate + 1.96*std.error)


## Create plot and export tiff to JOP specifications. 
plot_df <- 
  est_df %>%
  filter(!(model %in% c("Model 4", "Model 3", "Model 2") & 
             term == "Black Proposer")) %>%
  mutate(group = ifelse(model == "Model 4" & term == "Black Proposer x EPI",
                        "Triple-diff for prejudiced whites", term),
         group = ifelse(model == "Model 4" & term == "Black Proposer x RRI",
                        "Triple-diff for resentful whites", group),
         group = ifelse(model == "Model 3" & term == "Black Proposer x EPI",
                        "Diff-in-diff for prejudiced whites", group),
         group = ifelse(model == "Model 2" & term == "Black Proposer x RRI",
                        "Diff-in-diff for resentful whites", group),
         group = ifelse(model == "Model 1", "Main effect for all whites", 
                        group),
         group = factor(group, 
                        levels = rev(c("Main effect for all whites",
                                   "Diff-in-diff for resentful whites",
                                   "Diff-in-diff for prejudiced whites",
                                   "Triple-diff for resentful whites",
                                   "Triple-diff for prejudiced whites")),
                        labels = rev(c("M1.1: Main effect \n for all whites",
                                   "M1.2: Diff-in-diff for \n resentful whites",
                                   "M1.3: Diff-in-diff for \n prejudiced whites",
                                   "M1.4: Diff-in-diff for \n resentful whites",
                                   "M1.4: Diff-in-diff for \n prejudiced whites"))))

## Print results
plot_df %>% 
  dplyr::select(term, estimate, std.error, statistic, p.value, model) %>% 
  mutate_if(is.numeric, round, 2) 


plot_coef_ug <- 
  plot_df %>%
  ggplot(., aes(x = estimate, y = group)) + 
  geom_vline(xintercept = 0, size = 0.4, lty = 1) + 
  geom_errorbarh(aes(xmin = li95, xmax = ui95), height = 0,
                 size = 0.65, color = "black", alpha = 0.8) + 
  geom_errorbarh(aes(xmin = li90, xmax = ui90), 
                 size = 1.25, color = "black", height = 0) + 
  geom_point(size = 2, col = "black", fill = "black", show.legend = FALSE,
             pch = 22) +
  scale_x_continuous(limits = c(-0.05, 0.05),
                     labels = scales::percent_format(accuracy = 1)) +
  theme_kyle() + 
  xlab("Change in likelihood of offer acceptance") +
  ylab("") +
  theme(
    axis.line.x = element_line(size = 0.5, colour = "black"),
    axis.line.y = element_line(size = 0.5, colour = "black"),
    legend.position = "none",
    legend.key.size = unit(1, "line"),
    legend.title = element_text(size = 10, family = "serif"),
    legend.text = element_text(size = 10, family = "serif"),
    legend.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "cm"),
    plot.margin = unit(c(t = 0.2, r = 0.2, b = 1, l = 0.2), "lines"),
    strip.text.x = element_text(size = 12, family = "serif"),
    axis.title = element_text(size = 12, family = "serif"),
    axis.text.x = element_text(size = 12, family = "serif"),
    axis.text.y = element_text(size = 12, family = "serif"),
    axis.title.y = element_blank(),
    axis.ticks.x = element_line(colour = "black"),
    axis.ticks.y = element_blank()
  )

ggsave(paste0("output/fig1.tiff"), plot_coef_ug, height = 4, width = 5, 
       units = "in", dpi = 800)

#####-------------Study 2 Fig. 2 -----------------------------------------#####
fit1_eval <- felm(offer_fairness ~ round_racepair | factor(round) + 
                    as.factor(offer_amount) | 0 | pid, data = evals_ssi)

fit2_eval <- felm(offer_fairness ~ round_racepair*symbolic_binary | factor(round) + 
                    as.factor(offer_amount) | 0 | pid, data = evals_ssi)

fit3_eval <- felm(offer_fairness ~ round_racepair*overt_binary | factor(round) + 
                    as.factor(offer_amount) | 0 | pid, data = evals_ssi)

fit4_eval <- felm(offer_fairness ~ round_racepair*overt_binary + 
                    round_racepair*symbolic_binary | factor(round) + 
                    as.factor(offer_amount) | 0 | pid, data = evals_ssi)


# Compute BW-WB differences. BW-WB < 0 implies blacks perceived as less fair to 
# whites than whites are to blacks
base_est <- 
  summary(
    glht(fit1_eval, 
         linfct = c("round_racepairBW - round_racepairWB = 0")
         )
    )

rr_est <- 
  summary(
    glht(fit2_eval, 
         linfct = c("round_racepairBW:symbolic_binary - round_racepairWB:symbolic_binary = 0")
         )
    )

ep_est <- 
  summary(
    glht(fit3_eval, 
         linfct = c("round_racepairBW:overt_binary - round_racepairWB:overt_binary = 0")
         )
    )


rr_est2 <- 
  summary(
    glht(fit4_eval, 
         linfct = c("round_racepairBW:symbolic_binary - round_racepairWB:symbolic_binary = 0")
       )
    )

ep_est2 <- 
  summary(
    glht(fit4_eval, 
         linfct = c("round_racepairBW:overt_binary - round_racepairWB:overt_binary = 0")
       )
    )

diff_df <- 
  rbind(as.numeric(c(coef(base_est), sqrt(vcov(base_est)), 
                     base_est$test$tstat, base_est$test$pvalues)),
        as.numeric(c(coef(rr_est), sqrt(vcov(rr_est)),
                     rr_est$test$tstat, rr_est$test$pvalues)),
        as.numeric(c(coef(ep_est), sqrt(vcov(ep_est)),
                     ep_est$test$tstat, ep_est$test$pvalues)),
        as.numeric(c(coef(rr_est2), sqrt(vcov(rr_est2)),
                     rr_est2$test$tstat, rr_est2$test$pvalues)),
        as.numeric(c(coef(ep_est2), sqrt(vcov(ep_est2)),
                     ep_est2$test$tstat, ep_est2$test$pvalues))) %>% 
  data.frame() %>% 
  rename(estimate = X1, std.error = X2, statistic = X3, p.value = X4) %>% 
  mutate(model = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 4"),
         term = c("round_racepairBWWB", 
                  "round_racepairBWWB:symbolic_binary",
                  "round_racepairBWWB:overt_binary", 
                  "round_racepairBWWB:symbolic_binary",
                  "round_racepairBWWB:overt_binary"))


est_df <- 
  bind_rows(tidy(fit1_eval) %>% mutate(model = "Model 1"),
            tidy(fit2_eval) %>% mutate(model = "Model 2"),
            tidy(fit3_eval) %>% mutate(model = "Model 3"),
            tidy(fit4_eval) %>% mutate(model = "Model 4"),
            diff_df) %>% 
  filter(!(term %in% c("overt_binary", "symbolic_binary"))) %>% 
  mutate(term = factor(term, levels = c("round_racepairBB",
                                        "round_racepairBW",
                                        "round_racepairWB",
                                        "round_racepairBWWB",
                                        "round_racepairBB:symbolic_binary",
                                        "round_racepairBW:symbolic_binary",
                                        "round_racepairWB:symbolic_binary",
                                        "round_racepairBWWB:symbolic_binary",
                                        "round_racepairBB:overt_binary",
                                        "round_racepairBW:overt_binary",
                                        "round_racepairWB:overt_binary",
                                        "round_racepairBWWB:overt_binary"),
                       labels = c("BB - WW", 
                                  "BW - WW", 
                                  "WB - WW",
                                  "BW - WB",
                                  "(BB - WW) x RRI", 
                                  "(BW - WW) x RRI", 
                                  "(WB - WW) x RRI",
                                  "(BW - WB) x RRI",
                                  "(BB - WW) x EPI", 
                                  "(BW - WW) x EPI", 
                                  "(WB - WW) x EPI",
                                  "(BW - WB) x EPI"))) %>% 
  mutate(li90 = estimate - 1.645*std.error,
         ui90 = estimate + 1.645*std.error,
         li95 = estimate - 1.96*std.error,
         ui95 = estimate + 1.96*std.error)



## Create plot and export tiff to JOP specifications. 
plot_df <- 
  est_df %>%
  filter(!(model %in% c("Model 4", "Model 3", "Model 2") & 
             term %in% c("BB - WW", "BW - WW", "WB - WW"))) %>% 
  mutate(group = ifelse(model == "Model 4" & term == "(BW - WW) x EPI",
                        "Triple-diff for prejudiced whites", term),
         group = ifelse(model == "Model 4" & term == "(WB - WW) x EPI",
                        "Triple-diff for prejudiced whites", group),
         group = ifelse(model == "Model 4" & term == "(BB - WW) x EPI",
                        "Triple-diff for prejudiced whites", group),
         group = ifelse(model == "Model 4" & term == "(BW - WB) x EPI",
                        "Triple-diff for prejudiced whites", group),
         group = ifelse(model == "Model 4" & term == "(BW - WW) x RRI",
                        "Triple-diff for resentful whites", group),
         group = ifelse(model == "Model 4" & term == "(WB - WW) x RRI",
                        "Triple-diff for resentful whites", group),
         group = ifelse(model == "Model 4" & term == "(BB - WW) x RRI",
                        "Triple-diff for resentful whites", group),
         group = ifelse(model == "Model 4" & term == "(BW - WB) x RRI",
                        "Triple-diff for resentful whites", group),
         group = ifelse(model == "Model 3" & term == "(BW - WW) x EPI",
                        "Diff-in-diff for prejudiced whites", group),
         group = ifelse(model == "Model 3" & term == "(WB - WW) x EPI",
                        "Diff-in-diff for prejudiced whites", group),
         group = ifelse(model == "Model 3" & term == "(BB - WW) x EPI",
                        "Diff-in-diff for prejudiced whites", group),
         group = ifelse(model == "Model 3" & term == "(BW - WB) x EPI",
                        "Diff-in-diff for prejudiced whites", group),
         group = ifelse(model == "Model 2" & term == "(BW - WW) x RRI",
                        "Diff-in-diff for resentful whites", group),
         group = ifelse(model == "Model 2" & term == "(WB - WW) x RRI",
                        "Diff-in-diff for resentful whites", group),
         group = ifelse(model == "Model 2" & term == "(BB - WW) x RRI",
                        "Diff-in-diff for resentful whites", group),
         group = ifelse(model == "Model 2" & term == "(BW - WB) x RRI",
                        "Diff-in-diff for resentful whites", group),
         group = ifelse(model == "Model 1" & term == "BW - WW", 
                        "Main effect for all whites", group),
         group = ifelse(model == "Model 1" & term == "WB - WW", 
                        "Main effect for all whites", group),
         group = ifelse(model == "Model 1" & term == "BB - WW", 
                        "Main effect for all whites", group),
         group = ifelse(model == "Model 1" & term == "BW - WB", 
                        "Main effect for all whites", group),
         group = factor(group, 
                        levels = rev(c("Main effect for all whites",
                                   "Diff-in-diff for resentful whites",
                                   "Diff-in-diff for prejudiced whites",
                                   "Triple-diff for resentful whites",
                                   "Triple-diff for prejudiced whites")),
                        labels = rev(c("M2.1: Main effect \n for all whites",
                                   "M2.2: Diff-in-diff for \n resentful whites",
                                   "M2.3: Diff-in-diff for \n prejudiced whites",
                                   "M2.4: Diff-in-diff for \n resentful whites",
                                   "M2.4: Diff-in-diff for \n prejudiced whites"))),
         trt = ifelse(term == "(WB - WW) x RRI", "WB-WW", term),
         trt = ifelse(term == "(BW - WW) x RRI", "BW-WW", trt),
         trt = ifelse(term == "(BB - WW) x RRI", "BB-WW", trt),
         trt = ifelse(term == "(BW - WB) x RRI", "BW-WB", trt),
         trt = ifelse(term == "(BW - WW) x EPI", "BW-WW", trt),
         trt = ifelse(term == "(WB - WW) x EPI", "WB-WW", trt),
         trt = ifelse(term == "(BB - WW) x EPI", "BB-WW", trt),
         trt = ifelse(term == "(BW - WB) x EPI", "BW-WB", trt),
         trt = ifelse(term == "BW - WW", "BW-WW", trt),
         trt = ifelse(term == "WB - WW", "WB-WW", trt),
         trt = ifelse(term == "BB - WW", "BB-WW", trt),
         trt = ifelse(term == "BW - WB", "BW-WB", trt),
         trt = factor(trt, levels = c("BW-WW", "WB-WW", "BW-WB", "BB-WW"),
                      labels = c("Black to White (BW-WW)",
                                 "White to Black (WB-WW)",
                                 "Intergroup Difference (BW-WB)",
                                 "Black Proposer, Black Responder (BB-WW)")))

## Print results
plot_df %>% 
  dplyr::select(term, estimate, std.error, statistic, p.value, model) %>% 
  mutate_if(is.numeric, round, 2) %>% 
  arrange(model)

evals_coefplot <- 
  plot_df %>%
  filter(trt != "Black Proposer, Black Responder (BB-WW)") %>% 
  ggplot(., aes(x = estimate, y = group)) + 
  geom_vline(xintercept = 0, size = 0.4, lty = 1) + 
  geom_errorbarh(aes(xmin = li95, xmax = ui95), height = 0,
                 size = 0.65, color = "black", alpha = 0.8) + 
  geom_errorbarh(aes(xmin = li90, xmax = ui90), 
                 size = 1.25, color = "black", height = 0) + 
  geom_point(size = 2, col = "black", fill = "black", show.legend = FALSE,
             pch = 22) + 
  scale_x_continuous(limits = c(-3, 3)) +
  theme_kyle() + 
  facet_wrap(~trt) +
  xlab("Change in perceived fairness of proposed resource allocation") +
  ylab("") +
  theme(
    axis.line.x = element_line(size = 0.5, colour = "black"),
    axis.line.y = element_line(size = 0.5, colour = "black"),
    legend.position = "none",
    legend.key.size = unit(1, "line"),
    legend.title = element_text(size = 10, family = "serif"),
    legend.text = element_text(size = 10, family = "serif"),
    legend.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "cm"),
    plot.margin = unit(c(t = 0.2, r = 0.2, b = 1, l = 0.2), "lines"),
    strip.text.x = element_text(size = 10, family = "serif"),
    axis.title = element_text(size = 12, family = "serif"),
    axis.text.x = element_text(size = 12, family = "serif"),
    axis.text.y = element_text(size = 12, family = "serif"),
    axis.title.y = element_blank(),
    axis.ticks.x = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    panel.border = element_rect(colour = "black", size = 0.5)
  )



ggsave(paste0("output/fig2.tiff"), evals_coefplot, 
          height = 3.5, width = 8, dpi = 800, units = "in")

#####-------------Study 1 Descriptive Plot (Fig S1)-----------------------#####
ug_base_plot <- 
  mturk_ug %>% 
  mutate(Race = ifelse(proposer_race == "black", "Black proposer", 
                       "White proposer"),
         Race = factor(Race, levels = c("White proposer", 
                                        "Black proposer"))) %>% 
  group_by(offer_amount, proposer_race) %>%
  ggplot(., aes(x = offer_amount/100, y = accept_offer, group = Race, 
                shape = Race, colour = Race)) + 
  stat_smooth(se = FALSE, method = "loess") +
  scale_color_manual(guide = guide_legend(title = ""), 
                     values = c("black", "darkgrey"),
                     labels = c("White proposer", "Black proposer")) +
  theme_kyle() +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    y = "Likelihood of acceptance",
    x = "Percentage of endowment offered",
    caption = "Notes: Lines denote smoothed conditional means using local weighted regression (loess).") +
  guides(shape = FALSE, linetype = FALSE, group = FALSE) + 
  theme(legend.justification = c(1, 0), 
        legend.position = c(0.9, 0.35),
        text = element_text(size = 20, family = "serif"),
        axis.text = element_text(size = 18, family = "serif"),
        legend.title = element_blank(),
        axis.ticks.x = element_line(colour = "grey70", size = 0.1),
        axis.ticks.y = element_line(colour = "grey70", size = 0.1),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        panel.grid.major.y = element_line(colour = "grey70", size = 0.1))

ggplot2::ggsave(paste0("output/figS1.png"), ug_base_plot, height = 6, width = 8, 
                units = "in")

#####-------------Study 2 Descriptive Plot (Fig S2)-----------------------#####
fair_white_responder <- 
  evals_ssi %>% 
  filter(responder_race == "W", !is.na(offer_amount), !is.na(offer_fairness),
         offer_amount %in% c(10:40)) %>%
  mutate(responder_race = droplevels(responder_race)) %>% 
  ggplot(., aes(y = offer_fairness, x = factor(offer_amount), 
                color =  proposer_race, group = proposer_race, 
                linetype = proposer_race)) + 
  stat_smooth(method = "loess", se =  FALSE) + 
  scale_x_discrete(breaks =  seq(10, 40, 10), 
                   labels = paste0(seq(10, 40, 10), "%"))  +
  scale_linetype_manual(guide = guide_legend(title = ""), 
                        values = c("dashed", "solid"),
                        labels = c("Black Proposer", "White Proposer")) +
  scale_color_manual(guide = guide_legend(title = ""), 
                     values = c("darkgrey", "black"),
                     labels = c("Black Proposer", "White Proposer")) +
  ylab("Perceived Fairness of Offer") + 
  xlab("Percentage of endowment offered") +
  ggtitle("White Responder") +
  theme_kyle() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())


fair_black_responder <- 
  evals_ssi %>% 
  filter(responder_race == "B", !is.na(offer_amount), !is.na(offer_fairness),
         offer_amount %in% c(10:40)) %>%
  mutate(responder_race = droplevels(responder_race)) %>% 
  ggplot(., aes(y = offer_fairness, x = factor(offer_amount), 
                color =  proposer_race, group = proposer_race, 
                linetype = proposer_race)) + 
  stat_smooth(method = "loess", se =  FALSE) + 
  scale_x_discrete(breaks =  seq(10, 40, 10), 
                   labels = paste0(seq(10, 40, 10), "%"))  +
  scale_linetype_manual(guide = guide_legend(title = ""), 
                        values = c("dashed", "solid"),
                        labels = c("Black Proposer", "White Proposer")) +
  scale_color_manual(guide = guide_legend(title = ""), 
                     values = c("darkgrey", "black"),
                     labels = c("Black Proposer", "White Proposer")) +
  ylab("Perceived Fairness of Offer") + 
  xlab("Percentage of endowment offered") +
  ggtitle("Black Responder") +
  theme_kyle() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())

fair_white_proposer <- 
  evals_ssi %>% 
  filter(proposer_race == "W", !is.na(offer_amount), !is.na(offer_fairness),
         offer_amount %in% c(10:40)) %>%
  mutate(responder_race = droplevels(responder_race)) %>% 
  ggplot(., aes(y = offer_fairness, x = factor(offer_amount), 
                color =  responder_race, group = responder_race, 
                linetype = responder_race)) + 
  stat_smooth(method = "loess", se =  FALSE) + 
  scale_x_discrete(breaks =  seq(10, 40, 10), 
                   labels = paste0(seq(10, 40, 10), "%"))  +
  scale_linetype_manual(guide = guide_legend(title = ""), 
                        values = c("dashed", "solid"),
                        labels = c("Black Responder", "White Responder")) +
  scale_color_manual(guide = guide_legend(title = ""), 
                     values = c("darkgrey", "black"),
                     labels = c("Black Responder", "White Responder")) +
  ylab("Perceived Fairness of Offer") + 
  xlab("Percentage of endowment offered") +
  ggtitle("White Proposer") +
  theme_kyle() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())


fair_black_proposer <- 
  evals_ssi %>% 
  filter(proposer_race == "B", !is.na(offer_amount), !is.na(offer_fairness),
         offer_amount %in% c(10:40)) %>%
  mutate(responder_race = droplevels(responder_race)) %>% 
  ggplot(., aes(y = offer_fairness, x = factor(offer_amount), 
                color =  responder_race, group = responder_race, 
                linetype = responder_race)) + 
  stat_smooth(method = "loess", se =  FALSE) + 
  scale_x_discrete(breaks =  seq(10, 40, 10), 
                   labels = paste0(seq(10, 40, 10), "%"))  +
  scale_linetype_manual(guide = guide_legend(title = ""), 
                        values = c("dashed", "solid"),
                        labels = c("Black Responder", "White Responder")) +
  scale_color_manual(guide = guide_legend(title = ""), 
                     values = c("darkgrey", "black"),
                     labels = c("Black Responder", "White Responder")) +
  labs(y = "Perceived Fairness of Offer",
       x = "Percentage of endowment offered",
       title = "Black Proposer") +
  theme_kyle() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank()) 

p <- plot_grid(fair_white_responder, fair_black_responder, 
               fair_white_proposer, fair_black_proposer, 
               align = "v", nrow = 2) 

p <- 
  ggdraw(add_sub(p, label = "Notes: Lines denote smoothed conditional means using local weighted regression (loess)."))

save_plot(paste0("output/figS2.png"), p, base_height = 8, base_width = 10)


#####-------------Association Plot (Fig. S4)------------------------------#####
stacked_surveys <- 
  evals_ssi %>% 
  distinct(pid, .keep_all = TRUE) %>% 
  mutate(Study = "Study 2") %>% 
  bind_rows(mturk_ug %>% 
              distinct(responder_id, .keep_all = TRUE) %>% 
              mutate(Study = "Study 1"))

scatters <- 
  ggplot(stacked_surveys, aes(y = overt_index, x = symbolic_index, 
                              color = sample)) +
  geom_jitter(alpha = 0.35, shape = 19, size = 1.5) + 
  geom_smooth(se = FALSE, method = "loess") + 
  scale_color_manual(values = c("grey59", "black"),
                     labels = c("Study 1 (N = 741)", "Study 2 (N = 738)")) + 
  labs(x = "Racial Resentment Scale (0 to 1)", 
       y = "Explicit Prejudice Scale (-6 to 6)",
       caption = "Note: Panels A and B plot the distribution of the Explicit Prejudice and Resentment Scales.                         \n Panel C plots the bivariate relationships between Explicit Prejudice and Racial Resentment.              \n   Lines in panel C denote smoothed conditional means using local weighted regression (loess).           ") +
  scale_y_continuous(limits = c(-6,6), breaks = seq(-6,6,3)) +
  theme_kyle() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())

hist_overt <- 
  ggplot(stacked_surveys, aes(overt_index, ..density.. , fill = sample)) +
  geom_histogram(binwidth = 0.5, alpha = 0.6, position = "identity") + 
  scale_fill_manual(values = c("grey59", "black"),
                    labels = c("Study 1 (N = 741)", "Study 2 (N = 738)")) + 
  xlab("Explicit Prejudice Scale (-6 to 6)") + ylab("") +
  scale_x_continuous(limits = c(-6,6), breaks = seq(-6,6,3)) +
  theme_kyle() + guides(fill = FALSE) +
  theme(legend.title = element_blank(),
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())

hist_symbolic <- 
  ggplot(stacked_surveys, aes(x=symbolic_index, ..density.., fill=sample)) +
  geom_histogram(alpha = 0.6, binwidth = 0.1, position = "identity") +
  scale_fill_manual(values = c("grey59", "black"),
                    labels = c("Study 1 (N = 741)", "Study 2 (N = 738)")) + 
  xlab("Racial Resentment Scale (0 to 1)") + ylab("")  +
  theme_kyle() + guides(fill=FALSE) +
  theme(legend.title = element_blank(),
        strip.background = element_blank(),
        axis.ticks = element_line(colour = "grey70", size = 0.1),
        panel.grid.major = element_line(colour = "grey70", size = 0.1),
        panel.grid.minor = element_blank())

p <- plot_grid(plot_grid(hist_overt, hist_symbolic, align = "v", ncol = 2, 
                         labels = c('A', 'B')), scatters, ncol = 1, 
               rel_heights = c(1, 1.3),
               labels = c('','C'))

save_plot(paste0("output/figS4.png"), p, base_height = 9, base_width = 8)

